#include "AlignerAccuracy.h"

AlignerAccuracy::AlignerAccuracy() {
	correctlyMapped = 0;
	incorrectlyMapped = 0;
	exactlyMapped = 0;
	exactlyNotMapped = 0;
	notMapped = 0;
}

AlignerAccuracy::~AlignerAccuracy() {

}

int AlignerAccuracy::findAccuracy(string refFileName,
		string chromoSortedSamFileName, int editDist) {

	string line, chromosome, mapLocation, readSequence, temp, temp2,
			mappedSequence, moreHits = "", singleHit, complementarySeq;
	stringstream lineConv, hitsConv, singleHitConv;
	int index, flag = 0;
	int lengthOfRead, chromoNum, mappedLocation, currentChromoNum = -1;
	ifstream readFile;
	notMapped = 0;

	RefDB r;
	r.loadRefFile(refFileName);		// Loading the reference genome

	struct mapping real, mapped;// real denotes simulated read, mapped is found by mapper

	readFile.open(chromoSortedSamFileName.c_str());

	if (!readFile) {
		cout << "Error in opening " << chromoSortedSamFileName << endl;
		return -1;
	}
	/*
	 * Parsing sam file sorted by order of chromosome - sorted to complete mapping chromosome-wise
	 * Assuming file written in standard SAM format (verified with mrFAST output)
	 *
	 *
	 *
	 *
	 */
	while (getline(readFile, line)) {
		if (line[0] != '@') {
			//replace(line.begin(), line.end(), '_', ' ');

			lineConv.str("");
			lineConv.str(line);

			index = 1;
			flag = 0;
			while (lineConv >> temp) {

				switch (index) {

				case 3:
					if (temp == "*")
						flag = 1;
					mapped.chromosome = temp;
					index++;
					break;

				case 4:
					if (temp == "*")
						flag = 1;
					mapped.location = temp;
					index++;
					break;

				case 10:
					if (temp == "*")
						flag = 1;
					real.sequence = temp;
					index++;
					break;

				default:
					if (temp.size() > 2) {
						if (temp.substr(0, 3) == "XA:")
							moreHits = temp;
					}
					index++;
					break;
				}
			}

			/*
			 *  converting chromosome in string format to integer
			 *  Assuming header available in SAM output file
			 *  and chromosome in SAM header in same order
			 *  as in reference genome which is loaded above.
			 *
			 */

			if (flag != 1) {

				chromoNum = convertChromo(mapped.chromosome);
				//cout << mapped.chromosome << "::" << chromoNum << endl;

				// Checking if next chromosome is to be loaded
				if (currentChromoNum != chromoNum) {

					r.loadChromo(chromoNum);
					currentChromoNum = chromoNum;

				}

				lengthOfRead = real.sequence.size();

				stringstream(mapped.location) >> mappedLocation;

				mapped.sequence = r.getRefSeq(mappedLocation - 1,
						lengthOfRead + editDist);

				if (checkMapping(real.sequence, mapped.sequence, editDist))
					correctlyMapped++;
				else
					incorrectlyMapped++;

				//Checking for moreHits

				if (moreHits != "") {

					moreHits = moreHits.erase(0, 5);
					replace(moreHits.begin(), moreHits.end(), ';', ' ');

					hitsConv.str("");
					hitsConv.str(moreHits);
					while (hitsConv >> singleHit) {

						replace(singleHit.begin(), singleHit.end(), ',', ' ');
						singleHitConv.str("");
						singleHitConv.str(singleHit);

						singleHitConv >> temp2;
						mapped.chromosome = temp2;

						chromoNum = convertChromo(mapped.chromosome);
						//cout << mapped.chromosome << "::" << chromoNum << endl;

						// Checking if next chromosome is to be loaded
						if (currentChromoNum != chromoNum) {

							r.loadChromo(chromoNum);
							currentChromoNum = chromoNum;
						}

						singleHitConv >> mappedLocation;

						if (mappedLocation < 0) {

							mapped.sequence = r.getRefSeq(
									abs(mappedLocation) - 1,
									lengthOfRead + editDist);
							complementarySeq = getComplementary(
									mapped.sequence);

							if (checkMapping(real.sequence, complementarySeq,
									editDist))
								correctlyMapped++;
							else
								incorrectlyMapped++;
						} else {
							mapped.sequence = r.getRefSeq(
									abs(mappedLocation) - 1,
									lengthOfRead + editDist);
							if (checkMapping(real.sequence, mapped.sequence,
									editDist))
								correctlyMapped++;
							else
								incorrectlyMapped++;
						}

						singleHitConv.clear();
					}
					hitsConv.clear();
				}
				//no need of else, since do nothing is no moreHits
			} else
				notMapped++;

			lineConv.clear();

		}
	}
	readFile.close();
	return 0;
}

bool AlignerAccuracy::checkMapping(string realSequence, string mappedSequence,
		int editDistance) {
	cout << realSequence << endl << mappedSequence << endl;
	unsigned e = calEditDistance(mappedSequence, realSequence);
	cout << e << endl;
	if (e > editDistance) {
		cout << "unmapped" << endl;
		return false;
	}

	else {
		cout << "mapped" << endl;
		return true;
	}
}

string AlignerAccuracy::getComplementary(string mappedSeq) {

	char base;
	for (int i = 0; i < mappedSeq.size(); ++i) {
		base = mappedSeq[i];

		switch (base) {
		case 'A':
			mappedSeq[i] = 'T';
			break;
		case 'C':
			mappedSeq[i] = 'G';
			break;
		case 'T':
			mappedSeq[i] = 'A';
			break;
		case 'G':
			mappedSeq[i] = 'C';
			break;

		}
	}
	return mappedSeq;
}

unsigned int AlignerAccuracy::calEditDistance(string s1, string s2) {

	const size_t len1 = s1.size(), len2 = s2.size();
	vector<unsigned int> col(len1 + 1), prevCol(len1 + 1);
	unsigned int i, j, minVal;

	for (i = 0; i < prevCol.size(); i++)
		prevCol[i] = i;

	for (i = 0; i < len2; i++) {
		col[0] = i + 1;
		for (j = 0; j < len1; j++)
			col[j + 1] = min(min(1 + col[j], 1 + prevCol[1 + j]),
					prevCol[j] + (s2[i] == s1[j] ? 0 : 1));
		col.swap(prevCol);
	}
	minVal = prevCol[len1];
	for (i = len1; i > 0; i--) {
		if (prevCol[i - 1] < minVal)
			minVal = prevCol[i - 1];
	}
	return minVal;
}

unsigned int AlignerAccuracy::convertChromo(string chromosome) {
	return chromoMap[chromosome];
}

void AlignerAccuracy::buildChromoMap(string samFileName) {
	ifstream readFile;
	readFile.open(samFileName.c_str());
	string line, temp, chromo;
	stringstream lineConv;
	unsigned int index = 0;

	while (getline(readFile, line)) {

		if (line[0] != '@')
			break;

		lineConv.str("");
		lineConv.str(line);

		while (lineConv >> temp) {
			if ((temp.size() > 2) && (temp.substr(0, 3) == "SN:")) {
				chromo = temp.erase(0, 3);
				chromoMap[chromo] = index;
				index++;

				//cout << chromo << endl;
			}
		}
		lineConv.clear();
	}
	readFile.close();
}

void AlignerAccuracy::findCorrectness(string readSortedSamFileName,
		int errorThreshold) {
	cout << "inside findCorrectness" << endl;
	ifstream readFile;
	readFile.open(readSortedSamFileName.c_str());
	struct readSequence real, mapped, current;

	bool firstRun = true, matched = false;
	string line, mappedChromosome, mappedLocation, temp, trueChromosome,
			trueLocation, temp2, moreHits = "", singleHit;
	stringstream lineConv, readConv, hitsConv, singleHitConv;
	int index, idx, flag;
	unsigned long int lengthOfRead, chromoNum;
	notMapped = 0;

	while (getline(readFile, line)) {

		if (line[0] != '@') {
			cout << line << endl;
			//replace(line.begin(), line.end(), '_', ' ');

			lineConv.str("");
			lineConv.str(line);

			index = 1;
			flag = 0;
			while (lineConv >> temp) {
				switch (index) {

				case 1:
					replace(temp.begin(), temp.end(), '_', ' ');
					readConv.str("");
					readConv.str(temp);

					idx = 1;
					while (readConv >> temp2) {

						switch (idx) {

						case 1:
							if (temp2 == "*")
								flag = 1;
							real.chromosome = temp2;
							idx++;
							break;

						case 2:
							if (temp2 == "*")
								flag = 1;
							real.location = temp2;
							idx++;
							break;

						default:
							idx++;
							break;
						}

					}
					readConv.clear();
					index++;
					break;

				case 3:
					if (temp == "*")
						flag = 1;
					mapped.chromosome = temp;
					index++;
					break;

				case 4:
					if (temp == "*")
						flag = 1;
					mapped.location = temp;
					index++;
					break;

				default:
					if (temp.size() > 2) {
						if (temp.substr(0, 3) == "XA:")
							moreHits = temp;
					}
					index++;
					break;

				}
			}
			if (flag != 1) {

				if (firstRun == true) {
					current.chromosome = real.chromosome;
					current.location = real.location;
					firstRun = false;
				}

				if (!((real.chromosome == current.chromosome)
						&& (real.location == current.location))) {
					// Transition to next read
					if (matched)
						exactlyMapped++;
					else
						exactlyNotMapped++;

					matched = false;

					current.chromosome = real.chromosome;
					current.location = real.location;
				}

				cout << real.chromosome << "::" << real.location << endl;
				cout << mapped.chromosome << "::" << mapped.location << endl;
				if (checkMapping(mapped, real, errorThreshold))
					matched = true;

				//Doing for moreHits

				if (moreHits != "") {

					moreHits = moreHits.erase(0, 5);
					replace(moreHits.begin(), moreHits.end(), ';', ' ');

					hitsConv.str("");
					hitsConv.str(moreHits);
					while (hitsConv >> singleHit) {

						replace(singleHit.begin(), singleHit.end(), ',', ' ');

						singleHitConv.str("");
						singleHitConv.str(singleHit);

						singleHitConv >> temp2;
						mapped.chromosome = temp2;

						singleHitConv >> temp2;
						mapped.location = temp2;

						cout << mapped.chromosome << "::" << mapped.location
								<< endl;

						if (checkMapping(mapped, real, errorThreshold))
							matched = true;

						singleHitConv.clear();
					}
					hitsConv.clear();
				}

				lineConv.clear();

			} else
				notMapped++;
		}
	}
	if (matched)
		exactlyMapped++;
	else
		exactlyNotMapped++;

	readFile.close();

}

bool AlignerAccuracy::checkMapping(struct readSequence mapped,
		struct readSequence real, int errorThreshold) {

	unsigned int mappedLocation, realLocation;
	string mappedChromoParameter = mapped.chromosome;
	string mappedChromosome;
	stringstream chromoConv;
	chromoConv.str("");

	replace(mappedChromoParameter.begin(), mappedChromoParameter.end(), '_',
			' ');
	chromoConv.str(mappedChromoParameter);
	chromoConv >> mappedChromosome;

	if (mappedChromosome == real.chromosome) {

		stringstream(mapped.location) >> mappedLocation;
		stringstream(real.location) >> realLocation;

		if ((mappedLocation <= realLocation + errorThreshold)
				&& (mappedLocation >= realLocation - errorThreshold))
			return true;
		else
			return false;

	} else
		return false;
}

